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Abstract. - We present a molecular dynamics test of the Central Limit Theorem (CLT) in a 
paradigmatic long-range-interacting many-body classical Hamiltonian system, the HMF model. 
We calculate sums of velocities at equidistant times along deterministic trajectories for different 
sizes and energy densities. We show that, when the system is in a chaotic regime (specifically, 
at thermal equilibrium), ergodicity is essentially verified, and the Pdfs of the sums appear to 
be Gaussians, consistently with the standard CLT. When the system is, instead, only weakly 
chaotic (specifically, along longstanding metastable Quasi-Stationary States), nonergodicity (i.e., 
discrepant ensemble and time averages) is observed, and robust (/-Gaussian attractors emerge, 
consistently with recently proved generalizations of the CLT. 



Introduction. — During recent years there has been 
an increasing interest in generalizations of the Central 
Limit Theorem (CLT). This theorem - so called because of 
its central position in theory of probabilities - has ubiqui- 
tous and important applications in several fields. It essen- 
tially states that a (conveniently scaled) sum of n — ► oo in- 
dependent (or nearly independent) random variables with 
finite variance has a Gaussian distribution. Understand- 
ingly, this theorem is not applicable to those complex 
systems where long-range correlations are the rule, such 
as those addressed by noncxtensive statistical mechan- 
ics [1,2]. Therefore, several papers [3-10] have recently 
discussed extensions of the CLT and their correspond- 
ing attractors. In this paper, following [5,6], we present 
several numerical simulations for a long-range Hamilto- 
nian system, namely the Hamiltonian Mean Field (HMF) 
model. This model is a paradigmatic one for classical 
Hamiltonian systems with long-range interactions which 
has been intensively studied in the last decade (see, for 
example, [6,12-22], and references therein). In [5] it was 
shown that the probability density of rescaled sums of iter- 
ates of deterministic dynamical systems (e.g., the logistic 
map) at the edge of chaos (where the Lyapunov exponent 
vanishes) violates the CLT. Here we study rescaled sums of 
velocities considered along deterministic trajectories in the 



HMF model. It is well known that, in this model, a wide 
class of out-of-cquilibrium initial conditions induce a vio- 
lent relaxation followed by a metastable regime character- 
ized by nearly vanishing (strictly vanishing in the thermo- 
dynamic limit) Lyapunov exponents, and glassy dynam- 
ics [15-17]. We exhibit that correlations and nonergodicity 
created along these Quasi-Stationary States (QSS) can be 
so strong that, when summing the velocities calculated 
during the deterministic trajectories of single rotors at 
fixed intervals of time, the standard CLT is no longer ap- 
plicable. In fact, along the QSS, q-Gaussian Pdfs emerge 
as attractors instead of simple Gaussian Pdfs, consistently 
with the recently advanced g-generalized CLT [4,5,9], and 
ensemble averages arc different from time averages. 

Numerical simulations. — The HMF model de- 
scribes a system of N fully-coupled classical incrtial XY 

spins (rotors) s7= (cos 9i,sin Qi) , i — 1, iV,with uni- 
tary module and mass [12, 13]. These spins can also 
be thought as particles rotating on the unit circle. The 
Hamiltonian is given by 
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where 6i (0 < 6>, < 2ir) is the angle and pi the conjugate 
variable representing the rotational velocity of spin i. 
The equilibrium solution of the model in the canonical 
ensemble predicts a second order phase transition from a 
high temperature paramagnetic phase to a low tempera- 
ture ferromagnetic one [12]. The critical temperature is 
T c = 0.5 and corresponds to a critical energy per parti- 
cle U c = E c /N — 0.75. The order parameter of this phase 
transition is the modulus of the average magnetization per 

spin defined as: M = (1/N) \ ^2f =1 si | ■ Above T c , the 
spins point in different directions and M ~ 0. Below T c , 
most spins are aligned (the rotators are trapped in a single 
cluster) and M ^ 0. The out-of equilibrium dynamics of 
the model is also very interesting. In a range of energy 
densities between U <E [0.5, 0.75], special initial conditions 
called water-bag (characterized by initial magnetization 
Mq = 1 and uniform distribution of the momenta) drive 
the system, after a violent relaxation, towards mctastablc 
QSS. The latter slowly decay towards equilibrium with 
a lifetime which diverges like a power of the system size 
N [14-16]. 

In this section we simulate the dynamical evolution of 
several HMF systems with different sizes and at different 
energy densities, in order to explore their behavior either 
inside or outside the QSS regime. For each of them, fol- 
lowing the prescription of the CLT, we construct proba- 
bility density functions of quantities expressed as a finite 
sum of stochastic variables. But in this case, following the 
procedure adopted in ref. [5] for the logistic map, we will 
select these variables along the deterministics time evolu- 
tions of the N rotors. More formally, we study the Pdf of 
the quantity y defined as 

1 n 

Vi = ~r l>j(0~ < Pi >) f° r 3 = 1, 2, -,N , (2) 

vn i=i 

where Pj(i), with i = l,2,...,n, are the velocities of the 
jth-iotor taken at fixed intervals of time 5 along the same 
trajectory. The latter are obtained integrating the HMF 
equations of motions (see [15] for details about these equa- 
tions and the integration algorithm adopted ). The quan- 
tity < pj >= (1/n) Y^i=\Pj(f) 1S the average of the pj(i)'s 
over the single trajectory. The product Sx n gives the total 
simulation time. Note that the variables y's are propor- 
tional to the time average of the velocities along the single 
rotor trajectories. In the following we will distinguish this 
kind of average, i.e. time average, from the standard en- 
semble average, where the average of the velocities of the 
TV rotators is calculated at a given fixed time and over 
many different realizations of the dynamics. The latter 
can also be obtained from eq.(2) considering the y's vari- 
ables with n = 1 and < p 3 >= 0. In general, although 
the standard CLT predicts a Gaussian shape for sum of 
n independent stochastic values strictly when n — ► oo, in 
practice a finite sum converges quite soon to the Gaussian 
shape and this, in absence of correlations, is certainly true 
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Fig. 1: Temperature time evolution for the HMF system, with 
N=100 and Ml initial conditions, for U = 0.69 and for U = 0.4. 
The presence of a QSS regime is visible only in the U = 0.69 
case, although a transient regime exist also for U = 0.4. 



at least for the central part of the distribution [24] . Typi- 
cally we will use in this section a sum of n = 50 values of 
velocities along the deterministic trajectories for each of 
the N rotors of the HMF system, though larger values of 
n were also considered. 

In the following we will show that, if correlations among 
velocities are strong enough and the system is weakly 
chaotic, CLT predictions are not verified and, consistently 
with recent generalizations of the CLT, <j-Gaussians ap- 
pear [3-5]. The latter are a generalization of Gaussians 
which emerge in the context of noncxtensive statistical 
mechanics [1,2] and are defined as 

G q (x) ^ A(l - (1 - q)px 2 ) 1/1 - q , (3) 

being q the so-called entropic index (for q = 1 one re- 
covers the usual Gaussian) , (3 another suitable parameter 
(characterizing the width of the distribution), and A a 
normalization constant (see also ref. [10] for a simple and 
general way to generate them) . In particular we will show 
in this section that: 

(i) at equilibrium, when correlations are weak and the 
system is strongly chaotic (hence ergodic) standard CLT is 
verified, and time average coincides with ensemble average 
(both corresponding Pdfs arc Gaussians, either in the limit 
n — > oo or 5 — > oo); 

(ii) in the QSS regime, where velocities are strongly cor- 
related and the system is weakly chaotic and nonergodic, 
the standard CLT is no longer applicable, and q-Gaussian 
attractors replace the Gaussian ones; in this regime en- 
semble averages do not agree with time averages. 

For all the present simulations, water-bag initial condi- 
tions with initial magnetizazion Mq = 1, usually referred 
as Ml, will be used. In general, several different realiza- 
tions of the initial conditions will be performed also for 
the time average Pdfs case, but only in order to have a 
good statistics for small values of TV (for N=50000, on the 
contrary, only one realization has been used: see fig. 7(b)). 
Finally, to allow a correct comparison with standard Gaus- 
sians (represented as dashed lines in all the figures) and 
q-Gaussians (represented as full lines), the Pdf curves were 



Title 




5 -4 -3 -2 -1 (.1 1 2 

(y -<y>)te 




TIME AVERAGE 


(0)t 


-N=100-U=0.4-NO QSS 




E — q-Gauss 


an(q=1.35) 1 


n=50 jf \ 




5=200 fl 









TIME AVERAGE (d) 

N=100-U=0.4-NOQSS 

q-Gaussian (q=1 .1 ) 



4 6 8 10 -10 




Fig. 2: Numerical simulations for the HMF model with N=100, 
U=0.4 and Ml initial conditions. No QSS are present for this 
energy value, (a) We plot here the Pdf of the single rotor 
velocities at the time t=40000 (ensemble average over 1000 
realizations), i.e. we plot the (normalized) variable y defined as 
in eq.(2) with n = 1. The shape is Gaussian since the system 
is at equilibrium. In the other figures (b), (c) and (d) we 
plot the time average Pdfs for the normalized variable y, with 
n = 50 but with different time intervals (5=100,200 and 1000), 
calculated over an increasing simulation time after a transient 
of 40000 time units. An average over 1000 different realizations 
of the initial conditions was also considered in order to have 
a good statistics. Even if we are at equilibrium, it is evident 
a strong dependence of the entropic index q of the (/-Gaussian 
fitting curve on the time interval 8 adopted. Anyway, a time 
interval 8 = 1000 is already sufficient to obtain a Gaussian- 
shaped Pdf. See text for further details. 



always normalized to unit area and unit variance, by sub- 
tracting from the ?/'s their average < y > and dividing by 
the correspondent standard deviation a (hence, the tradi- 
tional yjn scaling adopted in Eq. (2) is in fact irrelevant). 



The case N=100. We start the discussion of the nu- 
merical simulations for the HMF model considering a size 
N = 100 and two different energy densities, U = 0.4 and 
U = 0.69. In the first case no QSS exist, while in the 
second case QSS characterize the out-of-equilibrium dy- 
namics and correlations formed during the first part of the 
dynamics decay slowly while the system relaxes towards 
equilibrium [15,16]. With N = 100 this relaxation takes 
however a reasonable amount of time steps, thus one can 
easily study also the equilibrium regime. The situation 
is illustrated in fig. 1, where we show the time evolu- 
tion of the temperature - calculated as twice the average 
kinetic energy per particle - for the two energy densities 
considered, starting from Mq = 1 initial conditions. As ex- 
pected QSS are clearly visible only in the case U = 0.69, 
although a small transient regime exists also for the case 
U = 0.4 [15]. 
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Fig. 3: Numerical simulations for the HMF model, with N — 
100 and U — 0.69 and Ml initial conditions. We are in this case 
inside the QSS regime, (a) We plot here the Pdf of the single 
rotor velocities at time t = 100 (ensemble average over 1000 
realizations). The shape is not Gaussian, (b) Time average 
Pdf for the normalized variable y with n = 50 and with a 
time interval 8 = 40, calculated after a transient time of 100 
time units. An average over 1000 different realizations of the 
initial conditions was also considered in order to have a good 
statistics. The resulting shape is very different from that one 
shown in (a) and can be very well fitted with a g-Gaussian. 



N=100 and U=0.4- Here we discuss numerical simu- 
lations for the HMF model with size N = 100 and U = 0.4. 
In this case it has been shown in the past that the equi- 
librium regime is reached quite fast and is characterized 
by a very chaotic dynamics [12,13]. 

In fig. 2 a transient time of 40000 units has been per- 
formed before the calculations, so that the equilibrium is 
fully reached (see fig.l). In (a) we consider the ensemble 
average of the velocities, i.e. the y variables defined as in 
(2) with n = 1, at t = 40000 and taking 1000 different re- 
alizations of the initial conditions (events). The Pdf com- 
pares very well with the Gaussian curve (dashed line), as 
expected at equilibrium. On the other hand, we consider 
in (b), (c) and (d) the Pdfs for the variable y with n = 50 
and with different time intervals 5 over an increasing sim- 
ulation time at equilibrium. As previously explained, this 
procedure corresponds to performing a time average along 
the trajectory for all the rotors of the system. In this 
case only the central part of the curve exhibits a Gaussian 
shape. On the other hand, Pdfs have long fat tails which 
can be very well reproduced with q-Gaussians (full lines) . 
If one increases the time interval S going from 5 = 100 
(b), to 8 = 200 (c) and finally to 6 = 1000 (d), the tails 
tend to disappear, the entropic index q of the q-Gaussians 
decreases from q = 1.45 ± 0.05 towards q = 1 and the 
Pdf tends to the standard Gaussian. This means that, as 
expected, summed velocities are less and less correlated 
as S increases (see also ref. [5]) and therefore the assump- 
tions of the CLT are satisfied as well as its prediction. 
Notice that n = 50 terms and a time interval <5 = 1000 are 
sufficiently large to reach a Gaussian-shaped Pdf. This 
situation reminds similar observations in the analysis of 
returns in financial markets [24], or in turbulence [25]. 

N=100 and U—0.69. Let us to consider now numeri- 
cal simulations for the HMF model with size N = 100 and 
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Fig. 4: Numerical simulations for the HMF model, with N = 
100, U — 0.69 and Ml initial conditions. We are here after the 
QSS regime, (a) Pdf of the single rotor velocities at t=40000 
(ensemble average over 1000 realizations). The shape is Gaus- 
sian since the system is at equilibrium. In the other figures we 
plot the time average Pdfs for the normalized variable y, calcu- 
lated after a transient time of 40000. An average over 1000 dif- 
ferent realizations of the initial conditions was also considered 
in order to have a good statistics. In figs.(b-d) we considered 
n = 50 but with different time intervals, more precisely (5=100 
(b), 5=1000 (c) and 5=2000 (d), over an increasing simulation 
time at equilibrium, (e)-(f) Pdfs obtained by keeping fixed the 
value S = 100 and increasing the number n of velocities in 
the sum for getting y. More precisely, n = 5000 in (e) and 
n = 50000 in (f). It is clear that, both for S — > oo and n — * do, 
the Pdfs shape tends to a Gaussian. 



U = 0.69. In this case a QSS regime exists, but its char- 
acteristic lifetime is quite short since the noise induced 
by the finite size drives the system towards equilibration 
rapidly. However strong correlations, created by the Ml 
initial conditions, exist and their decay is slower than in 
the case U = 0.4. In fig. 3 we show in (a) the Pdf of 
the velocities calculated at t = 100 (i.e. at the beginning 
of the QSS regime). An ensemble average over 1000 real- 
izations was considered. The Pdf shows a strange shape 
which remains constant in the QSS, as already observed 
in the past [14], and which differs from both the Gaussian 
and the g-Gaussian curves. On the other hand, we show in 
(b) the Pdf of the variable y with n = 50 and 6 = 40, i.e. 
calculated over a total of 2000 time steps after a transient 
of 100 units, in order to stay inside the QSS temperature 
plateaux (see fig.l). In this case the system is weakly 
chaotic and non ergodic [15,16] and the numerical Pdf is 
reproduced very well by a g-Gaussian with q = 1.65±0.05. 



Although in this case we have used differents initial con- 
ditions also for time averages, these results provide a first 
indication that ensemble and time averages are inequiva- 
lent in the QSS regime. Note that, due to the shortness 
of the QSS plateaux, for ./V = 100 it is not possible to use 
greater values of 5 or n in the numerical calculations of 
the y's. 

In fig. 4 we repeat the previous simulations for N = 100 
and U = 0.69, but adopting a transient time of 40000 
steps, in order to study the behavior of the system af- 
ter the QSS regime. The ensemble average Pdf (over 
1000 realizations) of the single rotor velocities at the time 
t = 40000 is shown in (a) and indicates that equilibrium 
seems to have been reached. In fact the agreement with 
the standard Gaussian is almost perfect up to 10~ 4 . In the 
other figures we plot the time average Pdfs for the variable 
y with n = 50 and for different time intervals <5, as done for 
U = 0.4. More precisely 5=100 in (b), (5=1000 in (c) and 
5=2000 in (d). Again it is evident a strong dependence 
of the Pdf shapes on the time interval 6 adopted. In fact 
initially (b) the Pdf is well fitted by a g-Gaussian with a 
q = 1.65 ± 0.05, however increasing <5, in (c) and (d), the 
central part of the Pdf becomes Gaussian while tails are 
still present and can be well fitted by q-Gaussians with 
values of q that tend towards unity. However, at variance 
with the U = 0.4 case, in this case not even a time interval 
(5 = 2000 is sufficient to reach a complete Gaussian-shaped 
Pdf down to 10~ 4 : evidently the strong correlations char- 
acterizing the QSS regime decay very slowly even after it, 
making the equilibrium shown by the ensemble average 
Pdf in (a) only apparent. This means that full ergodicity, 
i.e., full equivalence between ensemble and time averages, 
is reached, in this case, only asymptotically. 

The last statements are confirmed by panels (e) and (f) 
of fig. 4, where the effect of increasing the number n of 
summed velocities, keeping fixed the value of 5, has been 
investigated. More precisely 3=100 and n = 5000 in (e) 
and n = 50000 in (f). As expected, the increment of n 
makes the Pdf closer to the Gaussian, essentially because 
the total time over which the sum is considered increases 
(for n = 50000 we cover a simulation time of 5 x 10 6 ) 
and therefore correlations become asymptotically weaker 
and weaker, thus finally satisfying the prediction of the 
standard CLT 

In order to study in more details the ensemble-time in- 
equivalence along the QSS regime in the next subsection 
we will increase the system size and discuss numerical re- 
sults for N = 5000 and N = 50000. 

N=5000 and N=50000 at U=0.69. In fig.5 we show 
the time evolution of the temperature for the cases N = 
5000 and N = 50000 at U = 0.69, always starting (as 
usual) from the Ml initial conditions. It is evident that, 
for both systems, the length of the QSS plateaux is very 
much greater than for N = 100. 

We discuss first numerical simulations done inside the 
QSS for N = 5000 and U = 0.69. 
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Fig. 5: Temperature time evolution for the HMF system, with 
U = 0.69, Ml initial conditions and for TV = 5000 and TV = 
50000. The presence of a long-lasting QSS regime is clearly 
visible in both the cases and the plateaux are very much larger 
then in the TV = 100 case. 



In fig. 6 we show in (a) the ensemble average Pdf of ve- 
locities calculated over 1000 realizations at t = 100, i.e. 
at the beginning of the QSS regime. Its shape, constant 
along the entire QSS, is clearly not Gaussian and looks 
similar to that of fig. 3 (a). In panels (b-d) we show the 
effect of increasing the number n of velocity terms in the 
y sum on the time average Pdfs, calculated using a fixed 
value of 5 = 100. An average over 200 different realiza- 
tions of the initial conditions was also considered in order 
to have good statistics. In this case only for n = 1000 a 
g-Gaussian, with q = 1.45 ± 0.05, emerges. This is most 
likely due not to the effective number of n used but, con- 
sistently with fig. 6, to the fact that when choosing a large 
n one is averaging over a larger interval of time and thus 
considers in a more appropriate way the average over the 
entire QSS regime. In any case the observed behavior 
goes in the opposite direction to the prescriptions of the 
standard CLT and to the trend shown in panels (e-f) of 
fig. 4. Indeed, increasing n, the Pdf tails do not vanish 
but become more and more evident, thus supporting even 
further the claim about the existence of a non-Gaussian at- 
tractor for the nonergodic QSS regime of the HMF model. 
Moreover, the results of fig. 6 confirm the robustness of the 
(/-Gaussian shape along the entire QSS plateaux and the 
inequivalence between ensemble and time averages in the 
metastable regime. 

Let us now definitively demonstrate this inequivalence 
considering the case N=50000 at U=0.69. In fig. 7 (a) we 
plot the ensemble average Pdf of the velocities calculated 
(over 100 different realizations) at t = 200, i.e. at the be- 
ginning of the QSS regime, and after a very long transient, 
at t = 250000 (full circles). In panel (b) we plot the time 
average Pdf for the normalized variable y with n = 5000 
and 6 = 100, after a transient of 200 time units and over a 
simulation time of 500000 units along the QSS. It is impor- 
tant to stress that in this case only one single realization 
of the initial conditions has been performed, realizing this 
way a pure time average. The shape of the time aver- 
age Pdf (b) results to be again a robust g-Gaussian, with 
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Fig. 6: Numerical simulations for the HMF model with TV = 
5000, U — 0.69 and Ml initial conditions, in the QSS regime. 
(a) Pdf of single rotor velocities at the time t = 100 (ensemble 
average over 1000 realizations), (b-d) Time average Pdfs for 
the normalized variable y, after a transient time of 100 units 
and considering increasing values of n with a fixed time interval 
8 = 100, i.e. considering an increasing total simulation time 
inside the QSS. An average over 200 different realizations of 
the initial conditions was also considered in order to have a 
good statistics. Only for n = 1000, i.e. when the entire QSS 
extension has been considered (see fig. 5), we get a very good 
(/-Gaussian shape. 



q = 1.4 ± 0.05, not only in the tails, but also in the center 
(see inset). The time average Pdf is completely different 
from the ensemble average Pdf of fig. 7(a) (that is also very 
robust over all the plateaux), thus confirming definitively 
the inequivalence between the two kind of averages and 
the existence of a g-Gaussian attractor in the QSS regime 
of the HMF model. These results indicate that standard 
statistical mechanics based on the ergodic hypothesis can- 
not be applied in this case, while a generalized version, 
like the g-statistics [1,2] is likely more suitable [16]. 

Conclusions. — The numerical simulations presented 
in this paper strongly indicate that dynamical correlations 
and ergodicity breaking, induced in the HMF model by the 
initial out-of equilibrium violent relaxation, are present 
along the entire QSS metastable regime and decay very 
slowly even after it. In particular, considering finite sums 
of n correlated variables (velocities in this case) selected 
with a constant time interval S along single rotor trajec- 
tories, allowed us to study this phenomenon in detail. In- 
deed, we numerically showed that, in the weakly chaotic 
QSS regime, (i) ensemble average and time average of ve- 
locities are inequivalent, hence the ergodic hypothesis is 
violated, (ii) the standard CLT is violated, and (iii) robust 
g-Gaussian attractors emerge. On the contrary, when no 
QSS exist, or at a very large time after equilibration, i.e., 
when the system is fully chaotic and ergodicity has been 
restored, the ensemble average of velocities results to be 
equivalent to the time average and one observes a conver- 
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Fig. 7: Numerical simulations for the HMF model for N=50000, 
U=0.69 and Ml initial conditions in the QSS regime, (a) Pdfs 
of single rotor velocities at the times t=200 and t=250000 (en- 
semble average over 100 realizations), (b) Time average Pdf 
for the variable y calculated over only one single realization in 
the QSS regime and after a transient time of 200 units. In 
this case we used 8 = 100 and n = 5000, in order to cover a 
very large portion of the QSS (see fig. 5). Again, a g-Gaussian 
reproduces very well the calculated Pdf both in the tails and 
in the central part (see inset). See text for further details. 



gence towards the standard Gaussian attractor. In this 
case, the predictions of CLT are satisfied, even if we have 
only considered a finite sum of stochastic variables. How 
fast this happens depends on the size N , on the number 
n of terms summed in the y variables and on the time 
interval 5 considered. 

These results are consistent with the recent q- 
generalized forms of the CLT discussed in the literature 
[3-6,9], and pose severe questions to the often adopted 
procedure of using ensemble averages instead of time av- 
erages. Nonergodicity in coupled many particle systems 
goes back to the famous FPU experiment [26], but in our 
case is due to the long-range nature of the interaction. 
More recently, nonergodicity was found in deterministic 
iterative systems exibiting subdiffusion [11], but also in 
real experiments of shear flows, with results that were fit- 
ted with Lorentzians, i.e., g-Gaussians with q — 2 [23]. 
The whole scenario reminds that found for the leptokurtic 
returns Pdf in financial markets [24] , or in turbulence [25] , 
among many other systems, and could probably explain 
why y-Gaussians appear to be ubiquitous in complex sys- 
tems. Finally, we would like to add that, although it is 
certainly nontrivial to prove analytically whether the at- 
tractor in the nonergodic QSS regime of the HMF model 
precisely is a q-Gaussian or not (analytical results, as well 
as numerical dangers, have been recently illustrated in 
ref. [8] for various models), our numerical simulations un- 
ambiguously provide a very strong indication towards the 
existence of a robust g-Gaussian attractor in the case con- 
sidered. This opens new ways to the possible application 
of the y-gcncralized statistics in long-range Hamiltonian 
systems which will be explored in future papers. 
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